출산율 - 국가적 차원에서의 고찰

서론

유엔인구기금 2021년 세계 인구 현황 보고서 ‘내 몸은 나의 것’에 실린 통계표를 보면, 한국 여성 1명이 평생 낳을 것으로 예상되는 아이 수를 뜻하는 ‘합계 출산율’은 1.1명이다. 이렇게 한국의 저출산 문제와 인구 고령화가 심각해지면서 매년 저출산을 해결하기 위한 예산은 2006년부터 15년간 225조원을 들였지만 여전히 심각해져만 가고 있다. 또한 최근 2021년 1월 1일 기점으로 한국도 낙태법이 폐지되면서 낙태에 대한 처벌이 끝났다. 이번 한국의 낙태죄 폐지로 인해 출산율이 더 감소하는 것은 아닐까?하는 우려와 근심이 많이 제기되었다. 그렇다면 과연 어떤 요인이 출산율에 유의미하게 영향을 끼치는지, 우리나라의 저출산 해결 정책이 올바른 방향으로 가고 있는지 알아보고자 하였다.

가설

1. 결혼 나이가 낮을수록 출산율이 높을 것이다.
2. 교육 수준이 낮을수록 출산율이 높을 것이다.
   
3. 자살률이 낮을수록 출산율이 높을 것이다.
   
4. 육아비용이 낮을수록 출산율이 높을 것이다.
   
5. 65세 이상 인구 비율이 낮을수록 출산율이 높을 것이다.
   
6. 지니계수가 낮을수록 출산율이 높을 것이다.
   
7. 행복지수가 높을수록 출산율이 높을 것이다.
   
8. 낙태 합법인 국가는 출산율이 낮을 것이다.
   
9. GDP가 낮을수록 출산율이 낮을 것이다.
   
10. 평등지수가 높을수록 출산율이 높을 것이다.

본론

1 데이터 전처리

1.1 필요 라이브러리 불러오기

1.2 데이터 불러오기

우리는 전체 13곳에서 데이터를 모았다. 그 중 8개는 웹사이트에서 크롤링 하였고, 5개는 csv파일을 다운받아 사용하였다. 또한 수업 중 배운 R패키지를 통해 트위터 정보를 크롤링한 후, 저장한 뒤 사용하였다.

변수명	설명
sex ratio	여자 백명당 남자 몇명
gdp growth rate	연간 gdp 성장율
gdp per capita	국민 1인당 gdp
agriculture	농업 경제 % of GVA
industry	산업 경제 % of GVA
services and other activity	서비스 산업 % of GVA
agriculture employed	농업 종사자 비율
industry employed	산업 종사자 비율
service employed	서비스 종사자 비율
unemployed	실업자 비율
agriculture production	농업 생산지수
food production	식품 생산지수
exports	수출
imports	수입
trade balance	무역수지
urban population	도시 인구 비율
refugees	1000명당 난민이나 UNHCR에서 고려해야할 사람 수
infant mortality	유아 사망률
health expenditure	건강에 쓰는 비용 % of GDP
women politician	국회 여자 비율
internet using	거주지 100 당 인터넷 쓰는 사람 수
threatened species	멸종위기종수
CO2 emission	이산화탄소 예상 배출량
energy production	에너지 생산량

# 출산율
url <- "https://worldpopulationreview.com/country-rankings/birth-rate-by-country"
html_source <- GET(url)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F)
df_birth_rate = as.data.frame(tabs)
df_birth_rate <- df_birth_rate[,-3]
colnames(df_birth_rate) <- c("Country", "Birth rate")
df_birth_rate[,1] <- tolower(df_birth_rate[,1])
df_birth_rate$Country <- gsub(" ", "", df_birth_rate$Country)


# 결혼 나이 
url_first_marriage <- "https://en.wikipedia.org/wiki/List_of_countries_by_age_at_first_marriage"
html_source <- GET(url_first_marriage)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F)
df_first_marriage = data.frame(Country = character(), Men = numeric(), Women = numeric(), 
                               Average = numeric(), AgeGap = numeric(), AgeRatio = numeric(),
                               Year = numeric(), Source = character())
for (i in 1:5){
  table = as.data.frame(tabs[i])
  colnames(table) = table[1,]
  table = table[-1,]
  colnames(df_first_marriage) <- colnames(table)
  df_first_marriage = rbind(df_first_marriage, table)
}
df_first_marriage <- df_first_marriage[, c(1,2,3)]
colnames(df_first_marriage) <- c("Country", "Men marry", "Women marry")
df_first_marriage$Country <- tolower(df_first_marriage$Country)
df_first_marriage$Country <- gsub(" ", "", df_first_marriage$Country)


# 교육수준 
url <- "https://en.wikipedia.org/wiki/Education_Index"
html_source <- GET(url)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F)
df_study = as.data.frame(tabs[1])
colnames(df_study) = df_study[1,] ; df_study = df_study[-1,]
df_study <- df_study[, c(1, 31)]
colnames(df_study) <- c("Country", "education")
df_study$Country <- tolower(df_study$Country)
df_study$Country <- gsub(" ", "", df_study$Country)
df_study[91,1] <- "southkorea"
df_study$Country <- gsub("bosniaandherzegovina", "bosnia", df_study$Country)
df_study$Country <- gsub("unitedkingdom", "england", df_study$Country)
df_study$Country <- gsub("tanzania(unitedrepublicof)", "tanzania", df_study$Country)


# 자살률
url <- "https://en.wikipedia.org/wiki/List_of_countries_by_suicide_rate"
html_source <- GET(url)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F, header = F)
df_suicide = as.data.frame(tabs[1])
colnames(df_suicide) = df_suicide[1,] ; df_suicide = df_suicide[-1,]
df_suicide <- df_suicide[,-2]
colnames(df_suicide) <- c("Country", "Men suicide", "Women suicide")
df_suicide$Country <- tolower(df_suicide$Country)
df_suicide$Country <- gsub(" ", "", df_suicide$Country)
df_suicide$Country <- gsub("[[:punct:]]", "", df_suicide$Country)
df_suicide[22,1] <- "bosnia"
df_suicide[174,1] <- "england"


# 육아비용 
url <- "https://www.numbeo.com/cost-of-living/prices_by_country.jsp?displayCurrency=USD&itemId=224"
html_source <- GET(url)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F)
df_childcost = as.data.frame(tabs[2])
df_childcost = df_childcost[,-1]
colnames(df_childcost) <- c("Country", "child cost")
df_childcost$Country <- tolower(df_childcost$Country)
df_childcost$Country <- gsub(" ", "", df_childcost$Country)
df_childcost[45,1] <- "england"
df_childcost[61,1] <- "bosnia"


# 65세 이상 인구 비율
url <- "https://en.wikipedia.org/wiki/List_of_countries_by_age_structure"
html_source <- GET(url)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F)
df_oldratio = as.data.frame(tabs[1])
colnames(df_oldratio) = df_oldratio[1,] ; df_oldratio = df_oldratio[-1,]
df_oldratio = df_oldratio[-1,c(-2:-3)]
colnames(df_oldratio) <- c("Country", "old ratio")
df_oldratio$Country <- tolower(df_oldratio$Country)
df_oldratio$Country <- gsub(" ", "", df_oldratio$Country)
df_oldratio$Country <- gsub("bosniaandherzegovina", "bosnia", df_oldratio$Country)
df_oldratio$Country <- gsub("republicofthecongo", "congo", df_oldratio$Country)
df_oldratio$Country <- gsub("unitedkingdom", "england", df_oldratio$Country)


# 행복지수 
df_happiness <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/ahyoung/%EB%82%98%EB%9D%BC%EB%B3%84%20%ED%96%89%EB%B3%B5%EC%A7%80%EC%88%98.csv")
df_happiness <- df_happiness[,2:3]
colnames(df_happiness) <- c("Country", "happy score")
df_happiness$Country <- tolower(df_happiness$Country)
df_happiness$Country <- gsub(" ", "", df_happiness$Country)
df_happiness$Country <- gsub("bosniaandherzegovina", "bosnia", df_happiness$Country)
df_happiness$Country <- gsub("congo(brazzaville)", "congo", df_happiness$Country)
df_happiness$Country <- gsub("unitedkingdom", "england", df_happiness$Country)
df_happiness$Country <- gsub("southsudan", "sudan", df_happiness$Country)


# 지니계수 
df_gini <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/final_report/data_set/%EC%A7%80%EB%8B%88%EA%B3%84%EC%88%98.csv", fileEncoding = 'euc-kr')
df_gini[] <- lapply(df_gini, function(x){gsub("-", NA, x)})
rownames(df_gini) <- df_gini[,1]
df_gini[] <- lapply(df_gini, function(x){as.numeric(x)})
gini_mean <- round(rowMeans(df_gini, na.rm = T),2)
df_gini <- cbind(df_gini, gini_mean)
df_gini[,1] <- rownames(df_gini)
df_gini <- df_gini[,c(1,12)]
colnames(df_gini) <- c("Country", "gini_score")
df_gini$Country <- tolower(df_gini$Country)
df_gini$Country <- gsub(" ", "", df_gini$Country)


# 낙태 합법인 국가 
df_abortionlegal <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/abortionlegal.csv", header=FALSE)
df_abortionlegal <- df_abortionlegal[-219,-3]
df_abortionlegal <- as.data.frame(df_abortionlegal)
colnames(df_abortionlegal) = c("Country", "legal")
df_abortionlegal$legal[df_abortionlegal$legal != 'yes'] <- 'no'
df_abortionlegal$Country <- tolower(df_abortionlegal$Country)
df_abortionlegal$Country <- gsub(" ", "", df_abortionlegal$Country)


# 많은 데이터 
vary <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/final_report/data_set/%E1%84%82%E1%85%A1%E1%84%85%E1%85%A1%E1%84%83%E1%85%A6%E1%84%8B%E1%85%B5%E1%84%90%E1%85%A5.csv", header = FALSE)
colnames(vary) <- vary[1,]
vary <- vary[-1,]
df_vary <- as.data.frame(vary)
colnames(df_vary)[1] <- c("Country")
df_vary$Country <- tolower(df_vary$Country)
df_vary$Country <- gsub(" ", "", df_vary$Country)
df_vary$Country <- gsub("[[:punct:]]", "", df_vary$Country)
df_vary <- df_vary[,c(-2,-40,-41)]


# GDP
url <- "https://en.wikipedia.org/wiki/List_of_countries_by_GDP_(nominal)"
html_source <- GET(url)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F)
df_gdp = as.data.frame(tabs[3])
df_gdp <- df_gdp[5:194,-3]
colnames(df_gdp) <- c("Country", "gdp")
df_gdp$Country <- tolower(df_gdp$Country)
df_gdp$Country <- gsub(" ", "", df_gdp$Country)
df_gdp$Country <- gsub("[[:punct:]]", "", df_gdp$Country)
df_gdp[2,1] <- "china"
df_gdp$Country <- gsub("bosniaandherzegovina", "bosnia", df_gdp$Country)
df_gdp$Country <- gsub("solomonislands", "island", df_gdp$Country)
df_gdp$Country <- gsub("lebanon2020", "lebanon", df_gdp$Country)
df_gdp$Country <- gsub("moldovan8", "moldova", df_gdp$Country)
df_gdp$Country <- gsub("ukrainen5", "ukraine", df_gdp$Country)
df_gdp$Country <- gsub("unitedkingdom", "england", df_gdp$Country)
df_gdp$gdp<- gsub(",", "", df_gdp$gdp)


# 나라별 평등지수 
url <- "https://en.wikipedia.org/wiki/Global_Gender_Gap_Report"
html_source <- GET(url)
tabs <- readHTMLTable(rawToChar(html_source$content), stringsAsFactors=F)
df_equality = as.data.frame(tabs[3])
df_equality <- df_equality[-1,]
df_equality <- df_equality[-1,c(1,17)]
colnames(df_equality) = c("Country", "equality")
df_equality$Country <- tolower(df_equality$Country)
df_equality$Country <- gsub(" ", "", df_equality$Country)
df_equality$Country <- gsub("[[:punct:]]", "", df_equality$Country)
df <- full_join(df_gini, df_equality, by = 'Country')
df_equality$Country <- gsub("korearep", "southkorea", df_equality$Country)
df_equality$Country <- gsub("kyrgyzrepublic", "kyrgyzstan", df_equality$Country)
df_equality$Country <- gsub("bosniaandherzegovina", "bosnia", df_equality$Country)
df_equality$Country <- gsub("unitedkingdom", "england", df_equality$Country)
df_equality$Country <- gsub("democraticrepublicofthecongo", "congo", df_equality$Country)


# OECD 국가
df_oecd <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/final_report/data_set/OECD%20country.csv")
df_oecd <- cbind(df_oecd, 1)
colnames(df_oecd) <- c("Country", "OECD")

1.3 데이터 합치기, 이름바꾸기

우선 나라 이름을 기준으로 합쳐주었다. 한 나라에 대해서 다른 이름으로 된 나라들은 위에서 손수 바꿔주었다.

df <- inner_join(df_gini, df_abortionlegal, by = 'Country')
df <- inner_join(df, df_study, by = 'Country')
df <- inner_join(df, df_oldratio, by='Country')
df <- inner_join(df, df_happiness, by='Country')
df <- inner_join(df, df_childcost, by='Country')
df <- inner_join(df, df_suicide, by = "Country")
df <- inner_join(df, df_birth_rate, by = "Country")
df <- inner_join(df, df_gdp, by = 'Country')
df <- inner_join(df, df_equality, by='Country')
df <- inner_join(df, df_vary, by='Country')
df$Country <- str_to_title(df$Country)

colnames(df)[13:58] <- c("surface area", "population in thousands", "population density",
                         "sex ratio", "gdp_", "gdp growth rate", "gdp per capita", "agriculture",
                         "industry","services and other activity", "agriculture employed",
                         "industry employed","service employed", "unemployed", "labour force",
                         "agriculture production", "food production","exports", "imports",
                         "trade balance", "balance of paments", "population growth rate", "urban population",
                         "urban population growth rate", "fertility rate", "life expectancy at birth", 
                         "population age distribution", "international migrant", "refugees",
                         "infant mortality rate","health expenditure", "physicians",
                         "education government expenditure", "primary education", 
                         "secondary education", "tertiary education",
                         "women politician", "internet using", "threatened species",
                         "forested area", "CO2 emission", "energy production", "energy supply",
                         "water", "sanitation facilities", "net official")

1.4 데이터 변형

데이터들의 scale이 서로 다르기 때문에 모두 표준편차가1, 평균이 0인 상태로 만들어 주었다.
또한 본격적인 회귀분석에 앞서, 다중공선성 문제를 해결해 주기 위해 Birth rate를 종속변수로 하여 회귀분석을 한 뒤 다중공선성 문제가 있는 변수들을 삭제해 주었다.

# 낙태 합법 변수 재부호화
df$legal <- ifelse(df$legal == 'yes', 1, 0)

# 65 세 이상 변수 퍼센트 없애기
df$`old ratio` <- gsub("[[:punct:]]", "", df$`old ratio`)     # %없애주기
df$`old ratio` <- gsub("\u00A0", "", df$`old ratio`)          # 나머지 부분 없애주기

# 숫자형으로 변환
rownames(df) <- df$Country
df <- df[,-1]  # 나라이름 없애기
df[] <- lapply(df, function(x){as.numeric(x)})
no_na <- colSums(is.na(df))
na_col_list = vector()
for (i in 1:length(df)){
  if (no_na[[i]] != 0){
    na_col_list = append(na_col_list, i)
  }
}
df <- df[,na_col_list*-1]
df <- df[,-41]
df <- df[,c(-16, -32)]

# 표준화
df_std <- as.data.frame(scale(df))
head(sapply(df_std, sd))     # 표준편차

##  gini_score       legal   education   old ratio happy score  child cost 
##           1           1           1           1           1           1

head(sapply(df_std, mean))   # 평균

##                gini_score                     legal                 education 
## 0.00000000000000035060926 0.00000000000000007914230 0.00000000000000018397426 
##                 old ratio               happy score                child cost 
## 0.00000000000000001807053 0.00000000000000027642254 0.00000000000000005744489

# 다중공선성에 걸리는 변수 삭제
vif_model <- lm(`Birth rate` ~ ., data = df_std)
vif_value <- vif(vif_model) > 5
df_vif <- df_std[, vif_value]

2 국가별 출산율 분석

출산율은 연속적인 값이기 때문에 우리는 종속변수를 Birth rate로 한 뒤 다중선형 회귀 분석을 진행하였다.

2.1 다중선형 회귀분석

전진 선택법 사용하여 AIC가 낮아지는 쪽으로 변수를 선택하는 방법이다. 전진선택법을 통해 유의미한 변수들을 찾아보았다.

df_regression <- lm(`Birth rate` ~ ., data=df_vif)
df_regression_null <- lm(`Birth rate` ~ 1, data=df_vif)

step_df <- step(df_regression_null, scope=list(lower=df_regression_null, upper=df_regression),
                direction="forward")

summary(step_df) # 요약본 보기

## 
## Call:
## lm(formula = `Birth rate` ~ `old ratio` + `infant mortality rate` + 
##     `industry employed` + gini_score + `Women suicide`, data = df_vif)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.34898 -0.26033  0.00917  0.24655  1.42549 
## 
## Coefficients:
##                                         Estimate               Std. Error
## (Intercept)             -0.000000000000000009082  0.051711514739097380611
## `old ratio`             -0.491333196552296014303  0.083993501847600113286
## `infant mortality rate`  0.448512855986892722271  0.077955158697127033274
## `industry employed`     -0.144006310782061580644  0.056417257911671542248
## gini_score              -0.126978565252482467063  0.059534202151918375057
## `Women suicide`         -0.084507832070379720890  0.055480500329053147879
##                         t value    Pr(>|t|)    
## (Intercept)               0.000      1.0000    
## `old ratio`              -5.850 0.000000162 ***
## `infant mortality rate`   5.753 0.000000237 ***
## `industry employed`      -2.553      0.0130 *  
## gini_score               -2.133      0.0366 *  
## `Women suicide`          -1.523      0.1324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4418 on 67 degrees of freedom
## Multiple R-squared:  0.8183, Adjusted R-squared:  0.8048 
## F-statistic: 60.37 on 5 and 67 DF,  p-value: < 0.00000000000000022

round(lm.beta(step_df), 3)    # 계수확인 

##             `old ratio` `infant mortality rate`     `industry employed` 
##                  -0.491                   0.449                  -0.144 
##              gini_score         `Women suicide` 
##                  -0.127                  -0.085

모델을 보니 P값이 매우 낮으니, 유의미한 회귀 모델이다. AIC를 낮추는데 유의미한 변수는 5가지고, 회귀적으로 유의미한 변수는 4가지였다. 이 모델은 전체 데이터의 약 80%가량을 설명한다. 노인율과 산업고용율과 지니계수는 낮을 수록 출산율이 높았고, 유아 사망율은 높을 수록 출산율이 높았다. 유아사망율이 낮고, 산업고용율이 낮은 것은 개발도상국일 확률이 높다. 따라서 이들은 1차산업에 종사한 사람이 많기 때문에 노동력을 상승시키기 위해 더 많은 아를 낳는 것일 가능성이 있다. 따라서 인과관계는 없는것으로 생각된다. 또한 노인율 또한 인구구조 특성상 출산율이 높은 나라가 노인율이 작은 것이 대부분이다. 따라서 마찬가지로 인과관계는 없는 것으로 판단된다. 하지만 지니계수(소득불평등)가 낮을 수록 출산율이 높은 것은 빈부격차가 작을 수록 누구나 아이를 키우는데 힘들지 않기 때문에 출산율이 높은 것으로 볼 수 있다. 따라서 인과관계가 있을 수 있따. 더 자세한 결론은 밑에 결론부분에서 하겠다.

2.2 검증

회귀분석 결과를 신뢰할 수 있는지에 대한 검증을 시행하였다. 1번 그림은 다중공선성 문제를 살펴본 것이다. 결과는 문제가 없었다. 2번째는 qq곡선을 통해 정규성을 만족하는 지 여부를 본 것이다. 정규성을 만족한다고 볼 수 있다. 또한 세번째도 정규성을 검증하는 것이다. 마찮가지로 만족한다. 네 번째는 잔차들의 분포를 나타낸 것이다. 잔차들의 분포가 특정 모습을 띄지 않으므로, 본 회귀분석을 신뢰할 수 있다.

vif(step_df) # 다중공선성 확인

##             `old ratio` `infant mortality rate`     `industry employed` 
##                2.602115                2.241428                1.173976 
##              gini_score         `Women suicide` 
##                1.307279                1.135314

# sjplot
set_theme(axis.title.size = 1.0, axis.textsize = 1.0)
plot_model(step_df, type = "diag", wrap.labels=5)

## [[1]]

## 
## [[2]]

## `geom_smooth()` using formula 'y ~ x'

## 
## [[3]]

## 
## [[4]]

## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2,2))# plot 4개를 동시에 표시하기 위함
plot(step_df)

par(mfrow = c(1,1))

2.3 Elasticnet 분석

glmnet이라는 패키지를 이용하면, 단순히 정규화 뿐만 아니라 모델의 복잡도도 고려하여 회귀 모델을 만들 수 있다. 그래서 이번엔 elastic을 이용해 보겠다.

set.seed(1234) # 시드설정 
df_elastic <- cv.glmnet(as.matrix(df_vif[,-7]),
                        df_vif$`Birth rate`,
                        family = "gaussian", alpha = .5,
                        nfolds = 4, type.measure = "mse")

plot(df_elastic)      # 최적의 lambda값 확인

log(df_elastic$lambda.min)     # MSE가 가장 작을때의 lambda

## [1] -1.460013

log(df_elastic$lambda.1se)     # 가장 작을때 1표준편차 안에서 가장 간결한 모델일 때 lambda.

## [1] -0.9018106

plot(df_elastic$glmnet.fit, xvar = "lambda")   # 변수가 추가될 때마다 모수 추정값의 변화를 볼 수 있다.

coef.elastic <- coef(df_elastic, s = "lambda.min")[,1] 
coef.elastic       # 계수확인 

##                 (Intercept)                  gini_score 
## -0.000000000000000007108946  0.000000000000000000000000 
##                   education                   old ratio 
## -0.123567078742224403797856 -0.334688891616295614017673 
##                 happy score                 Men suicide 
##  0.000000000000000000000000  0.000000000000000000000000 
##               Women suicide                         gdp 
##  0.000000000000000000000000  0.000000000000000000000000 
##                surface area             gdp growth rate 
##  0.000000000000000000000000  0.000000000000000000000000 
##              gdp per capita                 agriculture 
##  0.000000000000000000000000  0.000000000000000000000000 
##                    industry services and other activity 
##  0.000000000000000000000000  0.000000000000000000000000 
##        agriculture employed           industry employed 
##  0.000000000000000000000000 -0.057301106232413190344488 
##                  unemployed      agriculture production 
##  0.000000000000000000000000  0.001137891967815737306532 
##             food production                     exports 
##  0.000000000000000000000000  0.000000000000000000000000 
##                     imports               trade balance 
##  0.000000000000000000000000  0.000000000000000000000000 
##                    refugees       infant mortality rate 
##  0.000000000000000000000000  0.301933059415346150000659 
##          threatened species                CO2 emission 
## -0.002939162053789053444253  0.000000000000000000000000 
##           energy production 
##  0.000000000000000000000000

mse.min.elastic <- df_elastic$cvm[df_elastic$lambda == df_elastic$lambda.min]
mse.min.elastic    # elestic 모델에서 최적의 mse값

## [1] 0.2763587

r2.min.elastic <- df_elastic$glmnet.fit$dev.ratio[df_elastic$lambda == df_elastic$lambda.min]
r2.min.elastic     # elestic 모델에서 최적의 R-squared값 

## [1] 0.7699512

# step_df 모델 수치
step_df_mse <- sum( (step_df$residuals)^2 )/nrow(df_vif)
step_df_mse

## [1] 0.1791634

결과적으로 보면 수치상 전진선택법의 모델이 mse값과 R^2값이 더 좋았지만, 모델의 복잡도 측면에서 본다면 위의 모델도 의미가 있는 모델이다.

2.4 Twitter 분석

#source(file = "TwitterAPIKey.R", echo = FALSE)

#childbirth <- enc2utf8("childbirth")
#data <- search_tweets(childbirth, n=10000, lang="en")
#data <- unique(data)
#data <- data[,c(1,2,3,4,5,6,7,67)]
#write_csv(data, "childbirthdata.csv")

data <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/final_report/data_set/childbirthdata.csv")

text <- str_replace_all(data$text, "@\\w+", "")
text <- str_replace_all(data$text, "[[:punct:]]", "") 
wordCorpus <- Corpus(VectorSource(text))
wordCorpus <- tm_map(wordCorpus, content_transformer(tolower))
wordCorpus <- tm_map(wordCorpus, removeWords, stopwords("english"))
wordCorpus <- tm_map(wordCorpus, removeWords, c("https"))
wordCorpus <- tm_map(wordCorpus, stripWhitespace)
wordCorpus <- tm_map(wordCorpus, removePunctuation)
wordCorpus <- tm_map(wordCorpus, removeNumbers)
wordCorpus <- tm_map(wordCorpus, removeWords, c("amp"))
wordCorpus <- tm_map(wordCorpus, removeWords, c("childbirth"))
wordCorpus <- tm_map(wordCorpus, removeWords, c("nhs", "realdonaldtrump", "tcoqzrnig"))
#wordCorpus[[1]]$content

pal <- brewer.pal(8, "Dark2")
wordcloud(words=wordCorpus, scale=c(4,0.5), max.words = 300,
          random.order = F, rot.per=0.35, use.r.layout = F, min.freq = 50, colors=pal)

encodesentiment <- function(x) {
  if(x <= -1){
    "1) very negative"
  }else if(x > -1 & x < 0){
    "2) negative"
  }else if(x > 0 & x < 1){
    "4) positive"
  }else if(x >= 1){
    "5) very positive"
  }else {
    "3) neutral"
  }
}

tweetsentiment <- get_sentiment(data$text, method= "syuzhet")
tweets <- cbind(data, tweetsentiment)
tweets$Sentiments <- sapply(tweets$tweetsentiment, encodesentiment)

qplot(tweets$tweetsentiment) + theme(legend.position = "none")+
  xlab("sentiment score") + ylab("number of tweets") +
  ggtitle("childbirth tweets by sentiment score")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(tweets, aes(Sentiments)) +
  geom_bar(fill="orange") +
  theme(legend.position = "none", axis.title = element_blank())+
  ylab("number of tweets") +
  ggtitle("tweets about childbirth")

#parenting <- enc2utf8("parenting")
#data2 <- search_tweets(parenting, n=10000, lang="en")
#data2 <- unique(data2)
#data2 <- data2[,c(1,2,3,4,5,6,7,67)]
#write_csv(data2, "parentingdata.csv")

data2 <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/final_report/data_set/parentingdata.csv")

text2 <- str_replace_all(data2$text, "@\\w+", "")
text2 <- str_replace_all(data2$text, "[[:punct:]]", "")
wordCorpus2 <- Corpus(VectorSource(text2))
wordCorpus2 <- tm_map(wordCorpus2, content_transformer(tolower))
wordCorpus2 <- tm_map(wordCorpus2, removeWords, stopwords("english"))
wordCorpus2 <- tm_map(wordCorpus2, removeWords, c("https"))
wordCorpus2 <- tm_map(wordCorpus2, stripWhitespace)
wordCorpus2 <- tm_map(wordCorpus2, removePunctuation)
wordCorpus2 <- tm_map(wordCorpus2, removeNumbers)
wordCorpus2 <- tm_map(wordCorpus2, removeWords, c("parenting"))
wordCorpus2 <- tm_map(wordCorpus2, removeWords, c("tcolukpux", "yoongi"))
wordCorpus2 <- tm_map(wordCorpus2, removeWords, c("httpstcolukpux"))
wordCorpus2 <- tm_map(wordCorpus2, removeWords, c("amp"))
#wordCorpus2[[1]]$content

set.seed(1234)
wordcloud(words=wordCorpus2, scale=c(4,0.5), max.words = 200,
          random.order = F, rot.per=0.35, use.r.layout = F, min.freq = 10, colors=pal)

tweetsentiment2 <- get_sentiment(data2$text, method= "syuzhet")
tweets2 <- cbind(data2, tweetsentiment2)
tweets2$Sentiments2 <- sapply(tweets2$tweetsentiment2, encodesentiment)

qplot(tweets2$tweetsentiment2) + theme(legend.position = "none")+
  xlab("sentiment score") + ylab("number of tweets") +
  ggtitle("parenting tweets by sentiment score")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(tweets2, aes(Sentiments2)) +
  geom_bar(fill="aquamarine4") +
  theme(legend.position = "none", axis.title = element_blank())+
  ylab("number of tweets") +
  ggtitle("tweets about parenting")

#childbirth_k <- enc2utf8("출산")
#Kdata <- search_tweets(childbirth_k, n=10000, lang="ko")
#Kdata <- unique(Kdata)
#Kdata <- Kdata[,1:13]
#write_csv(Kdata, "출산.csv")
Kdata <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/final_report/data_set/%E1%84%8E%E1%85%AE%E1%86%AF%E1%84%89%E1%85%A1%E1%86%AB.csv")

Kdata_text <- Kdata$text
Kdata_text <- sapply(Kdata_text, extractNoun, USE.NAMES = F)
Kdata_word <- unlist(Kdata_text)
Kdata_word <- Filter(function(x){nchar(x)>=2},Kdata_word)

Kdata_word <- gsub("\n", "", Kdata_word)
Kdata_word <- gsub("\r", "", Kdata_word)
Kdata_word <- gsub("https://", "", Kdata_word)
Kdata_word <- gsub("@", "", Kdata_word)
Kdata_word <- gsub("[[:punct:]]", "", Kdata_word)
Kdata_word <- str_replace_all(Kdata_word,"[A-Za-z0-9]","")
Kdata_word <- gsub("[ㄱ-ㅎ]", "", Kdata_word)
Kdata_word <- gsub("[ㅜ|ㅠ]", "", Kdata_word)
Kdata_word <- gsub("로", "", Kdata_word)
Kdata_word <- gsub("합사시켜", "", Kdata_word)
Kdata_word <- gsub("해서", "", Kdata_word)
Kdata_word <- gsub("출산", "", Kdata_word)
Kdata_word <- gsub("쿠오모가", "", Kdata_word)
Kdata_word <- gsub("까지", "", Kdata_word)
Kdata_word <- gsub("들이", "", Kdata_word)
Kdata_word <- gsub("그래서", "", Kdata_word)
Kdata_word <- gsub("그걸", "", Kdata_word)

head(Kdata_word)

## [1] "암수"   "마리"   "물리적" "연결"   "혈액"   "공유"

set.seed(1234)
wordcloud(words=Kdata_word, scale=c(5,0.5),
          min.freq = 30, random.order = F, max.words=200, family="AppleGothic", colors=pal)

senti_words_kr <- readr::read_delim("https://raw.githubusercontent.com/park1200656/KnuSentiLex/master/SentiWord_Dict.txt", delim='\t', col_names=c("term", "score"))

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   term = col_character(),
##   score = col_double()
## )

head(senti_words_kr)

## # A tibble: 6 x 2
##   term  score
##   <chr> <dbl>
## 1 (-;       1
## 2 (;_;)    -1
## 3 (^^)      1
## 4 (^-^)     1
## 5 (^^*      1
## 6 (^_^)     1

x <- duplicated(senti_words_kr$term)
senti_words_kr2 <- senti_words_kr[!x, ]
senti_dic_kr <- SentimentDictionaryWeighted(words = senti_words_kr2$term, 
                                            scores = senti_words_kr2$score)
senti_dic_kr <- SentimentDictionary(senti_words_kr2$term[senti_words_kr2$score > 0], 
                                    senti_words_kr2$term[senti_words_kr2$score < 0])

res_sentiment <- analyzeSentiment(Kdata_word,
                                  language="korean",
                                  rules=list("KoreanSentiment"=list(ruleSentiment, senti_dic_kr)),
                                  removeStopwords = F, stemming = F)
theme_set(theme_minimal(base_family = "AppleGothic"))

df <- data.frame(res_sentiment, Kdata_word)
df <- df[!(df$Kdata_word==""),]

df[is.na(df$KoreanSentiment),1] <- "neutral"
df[df$KoreanSentiment == 0,1] <- "neutral"
df[df$KoreanSentiment == 1,1] <- "positive"
df[df$KoreanSentiment == 0.125,1] <- "positive"
df[df$KoreanSentiment == 0.111111111111111,1] <- "positive"
df[df$KoreanSentiment == 0.0833333333333333,1] <- "positive"
df[df$KoreanSentiment < 0,1] <- "negative"

df_final <- df[!(df$KoreanSentiment=="neutral"),] #중립이 압도적으로 많아 빼고 시각화

ggplot(df_final, aes(x = KoreanSentiment)) + 
  geom_bar(stat = "count", width = 0.7, fill = "orange") + 
  theme_minimal() + ggtitle("tweets about childbirth in korea")

#childcare_k <- enc2utf8("육아")
#Kdata2 <- search_tweets(childcare_k, n=10000, lang="ko")
#Kdata2 <- unique(Kdata2)
#Kdata2 <- Kdata2[,1:13]
#write_csv(Kdata2, "육아.csv")

Kdata2 <- read.csv("https://raw.githubusercontent.com/jkworldchampion/Data_analytics/main/final_report/data_set/%E1%84%8B%E1%85%B2%E1%86%A8%E1%84%8B%E1%85%A1.csv")

Kdata_text2 <- Kdata2$text

Kdata_text2 <- sapply(Kdata_text2, extractNoun, USE.NAMES = F)
Kdata_word2 <- unlist(Kdata_text2)
Kdata_word2 <- Filter(function(x){nchar(x)>=2},Kdata_word2)

Kdata_word2 <- gsub("\n", "", Kdata_word2)
Kdata_word2 <- gsub("\r", "", Kdata_word2)
Kdata_word2 <- gsub("https://", "", Kdata_word2)
Kdata_word2 <- gsub("@", "", Kdata_word2)
Kdata_word2 <- gsub("[[:punct:]]", "", Kdata_word2)
Kdata_word2 <- str_replace_all(Kdata_word2,"[A-Za-z0-9]","")
Kdata_word2 <- gsub("[ㄱ-ㅎ]", "", Kdata_word2)
Kdata_word2 <- gsub("[ㅜ|ㅠ]", "", Kdata_word2)
Kdata_word2 <- gsub("로", "", Kdata_word2)
Kdata_word2 <- gsub("해서", "", Kdata_word2)
Kdata_word2 <- gsub("육아", "", Kdata_word2)
Kdata_word2 <- gsub("키선배님이", "", Kdata_word2)
Kdata_word2 <- gsub("놀토가", "", Kdata_word2)
Kdata_word2 <- gsub("호시는", "", Kdata_word2)
Kdata_word2 <- gsub("트니트니", "", Kdata_word2)
Kdata_word2 <- gsub("세븐", "", Kdata_word2)
Kdata_word2 <- gsub("그거", "", Kdata_word2)
Kdata_word2 <- gsub("갖다대니까", "", Kdata_word2)
Kdata_word2 <- gsub("하게", "", Kdata_word2)

set.seed(1234)
wordcloud(words=Kdata_word2, scale=c(4,0.5),
          min.freq = 40, random.order = F, max.words=200, family="AppleGothic", colors=pal)

res_sentiment2 <- analyzeSentiment(Kdata_word2,
                                   language="korean",
                                   rules=list("KoreanSentiment"=list(ruleSentiment, senti_dic_kr)),
                                   removeStopwords = F, stemming = F)

df2 <- data.frame(res_sentiment2, Kdata_word2)
df2 <- df2[!(df2$Kdata_word2==""),]

df2[is.na(df2$KoreanSentiment),1] <- "neutral"
df2[df2$KoreanSentiment == 0,1] <- "neutral"
df2[df2$KoreanSentiment == 1,1] <- "positive"
df2[df2$KoreanSentiment == 0.0555555555555556,1] <- "positive"
df2[df2$KoreanSentiment == 0.333333333333333,1] <- "positive"
df2[df2$KoreanSentiment == 0.142857142857143,1] <- "positive"
df2[df2$KoreanSentiment == 0.1,1] <- "positive"
df2[df2$KoreanSentiment == 0.0833333333333333,1] <- "positive"
df2[df2$KoreanSentiment == 0.0769230769230769,1] <- "positive"
df2[df2$KoreanSentiment == 0.0588235294117647,1] <- "positive"
df2[df2$KoreanSentiment == 0.05,1] <- "positive"
df2[df2$KoreanSentiment < 0,1] <- "negative"
df2[df2$KoreanSentiment == 0.0714285714285714,1] <- "negative"

df2_final <- df2[!(df2$KoreanSentiment=="neutral"),]

ggplot(df2_final, aes(x = KoreanSentiment)) + 
  geom_bar(stat = "count", width = 0.7, fill = "aquamarine4") + 
  theme_minimal() + ggtitle("tweets about parenting in korea")

3 OECD 회원국 분석

OECD는 특정한 가입 조건이 명시되어 있지는 않다. 하지만 OECD에 가입하려면 OECD 소속국 간의 허가를 받아야 속할 수 있다. 그래서 우리는 OECD 소속 국가만의 특징이 있을지 궁금해 졌다. 따라서 우리는 로지스틱 회귀분석으로 OECD국의 특징에 대해 알아보고자 한다.

3.1 로지스틱 회귀분석

OECD국가 일때 1이며 아닐때 0이다. 종속변수를 OECD국가로 한 뒤 로지스틱 회귀분석을 실시하였다.

Country <- rownames(df_vif)
df_vif <- cbind(df_vif, Country)
df_vif <- left_join(df_vif, df_oecd, by='Country')
df_vif$OECD[is.na(df_vif$OECD)] <- 0  # OECD여부에 따라 맞으면1, 아니면0
rownames(df_vif) <- df_vif[,28]
df_vif <- df_vif[,-28] # Country 삭제
df_vif$OECD <- as.numeric(df_vif$OECD)

# 다중공선성에 의한 변수 삭제
vif_model <- lm(`OECD` ~ ., data = df_vif)
vif_value <- vif(vif_model) > 10
df_vif_modify <- df_vif[, !vif_value]

df_log_reg <- glm(`OECD` ~ ., family = binomial, data = df_vif_modify, maxit = 100) # 변수 초과에 의해 정답률이 1이 나온다.

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

summary(df_log_reg)

df_log_reg_Null = glm(OECD ~ 1, family = binomial, data = df_vif_modify, maxit = 100)
summary(df_log_reg_Null)

# anova(df_log_reg_Null, df_log_reg, test="LRT")

step_df_vif_log_reg <- step(df_log_reg_Null, scope=list(lower=df_log_reg_Null, upper=df_log_reg),
                            direction="forward")

summary(step_df_vif_log_reg)

## 
## Call:
## glm(formula = OECD ~ `gdp per capita` + `old ratio` + refugees + 
##     `happy score` + `threatened species` + `industry employed` + 
##     unemployed, family = binomial, data = df_vif_modify, maxit = 100)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.54586  -0.10218  -0.00471   0.04423   1.85841  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)  
## (Intercept)           -1.0561     0.7486  -1.411   0.1583  
## `gdp per capita`       2.9799     1.6211   1.838   0.0660 .
## `old ratio`            1.6961     0.8029   2.113   0.0346 *
## refugees               2.5362     1.1460   2.213   0.0269 *
## `happy score`          3.9324     1.6149   2.435   0.0149 *
## `threatened species`   1.5092     0.7658   1.971   0.0487 *
## `industry employed`    1.9974     1.1339   1.762   0.0782 .
## unemployed             1.2554     0.9011   1.393   0.1636  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 99.536  on 72  degrees of freedom
## Residual deviance: 24.419  on 65  degrees of freedom
## AIC: 40.419
## 
## Number of Fisher Scoring iterations: 8

exp(coef(step_df_vif_log_reg))      # 계수를 살펴본다 

##          (Intercept)     `gdp per capita`          `old ratio` 
##            0.3478213           19.6852215            5.4528847 
##             refugees        `happy score` `threatened species` 
##           12.6318049           51.0281438            4.5233084 
##  `industry employed`           unemployed 
##            7.3697166            3.5091747

vif(step_df_vif_log_reg)            # 다중공선성 확인

##     `gdp per capita`          `old ratio`             refugees 
##             2.493158             2.011996             2.735693 
##        `happy score` `threatened species`  `industry employed` 
##             3.936967             1.932457             3.257611 
##           unemployed 
##             3.854877

로지스틱 결과 위에서 본 회귀계수는 오즈비이다(odds의 비율). 즉 OECD 국가가 더 선진국에 가까운 모습이 보인다.
GDP가 높고, 행복도가 높고, 산업고용율이 높을 수록 OECD일 확률이 높다.

3.1.1 모델 평가

df_log_reg_graph <- ROC(form = `OECD` ~ `old ratio`+
                          `infant mortality rate`+
                          `industry employed`+
                          `gini_score`+
                          `Women suicide`,
                        data = df_vif_modify,
                        plot = "ROC")

df_log_reg_graph$res[round(df_log_reg_graph$res$lr.eta,3) == 0.674,]

##                        sens      spec       pvp       pvn    lr.eta
## 0.673599504785467 0.8387097 0.9285714 0.1136364 0.1034483 0.6735995

df_log_reg_graph$AUC 

## [1] 0.9162826

df_log_reg_graph$lr

## 
## Call:  glm(formula = form, family = binomial, data = data)
## 
## Coefficients:
##             (Intercept)              `old ratio`  `infant mortality rate`  
##                -1.89486                  0.86574                 -4.21271  
##     `industry employed`               gini_score          `Women suicide`  
##                -0.02506                  0.96928                  0.07943  
## 
## Degrees of Freedom: 72 Total (i.e. Null);  67 Residual
## Null Deviance:       99.54 
## Residual Deviance: 54.49     AIC: 66.49

# caret
confusionMatrix(
  as.factor(ifelse(predict(step_df_vif_log_reg, type = "response")>0.674,1,0)),
  as.factor(df_log_reg$y),
  positive = '1')

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 41  4
##          1  1 27
##                                           
##                Accuracy : 0.9315          
##                  95% CI : (0.8474, 0.9774)
##     No Information Rate : 0.5753          
##     P-Value [Acc > NIR] : 0.00000000001086
##                                           
##                   Kappa : 0.858           
##                                           
##  Mcnemar's Test P-Value : 0.3711          
##                                           
##             Sensitivity : 0.8710          
##             Specificity : 0.9762          
##          Pos Pred Value : 0.9643          
##          Neg Pred Value : 0.9111          
##              Prevalence : 0.4247          
##          Detection Rate : 0.3699          
##    Detection Prevalence : 0.3836          
##       Balanced Accuracy : 0.9236          
##                                           
##        'Positive' Class : 1               
## 

F1_score <- 2*(0.9643*0.8710)/(0.9643+0.8710)
F1_score

## [1] 0.9152785

3.2 나이브 베이즈

# 그룹화 함수
mygroup = function (y,k=4){
  count = length(y)
  z = rank(y,ties.method = "min")
  return(floor((z-1)/(count/k))+1)
}
df_vif_factor <- data.frame()
df_vif_factor <- lapply(df_vif, function(x){mygroup(x,10)}) # 각각 10개씩 그룹으로 나눈다.
df_vif_factor <- as.data.frame(df_vif_factor)

df_vif_factor$OECD <- ifelse(df_vif_factor$OECD == '6', 1, 0)  # 종속도 변하기 때문에 다시 바꿔준다

# 데이터 분할
train_idx = sample(73, 73*3/4)      
df_train = df_vif_factor[train_idx,]
df_test = df_vif_factor[-train_idx,]
df_train_labels = df_vif_factor[train_idx,]$OECD
df_test_labels = df_vif_factor[-train_idx,]$OECD


df_train = df_train[,-length(df_train)] # OECD 열을 삭제

prop.table(table(df_train_labels))

## df_train_labels
##         0         1 
## 0.5925926 0.4074074

prop.table(table(df_test_labels))

## df_test_labels
##         0         1 
## 0.5263158 0.4736842

# 비슷하다.

# 모델 생성
df_classifier = naiveBayes(df_train, df_train_labels)

3.2.1 모델 평가

df_test_pred = predict(df_classifier, df_test)

CrossTable(df_test_pred, df_test_labels,
           prop.chisq = F, prop.c = F, prop.r = FALSE,
           dnn = c('predicted', 'actual'))

##    Cell Contents 
## |-------------------------|
## |                       N | 
## |         N / Table Total | 
## |-------------------------|
## 
## ==================================
##              actual
## predicted        0       1   Total
## ----------------------------------
## 0                9       4      13
##              0.474   0.211        
## ----------------------------------
## 1                1       5       6
##              0.053   0.263        
## ----------------------------------
## Total           10       9      19
## ==================================

confusionMatrix(as.factor(df_test_pred),
                as.factor(df_test_labels), 
                positive = '1')

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 0 1
##          0 9 4
##          1 1 5
##                                          
##                Accuracy : 0.7368         
##                  95% CI : (0.488, 0.9085)
##     No Information Rate : 0.5263         
##     P-Value [Acc > NIR] : 0.05192        
##                                          
##                   Kappa : 0.4633         
##                                          
##  Mcnemar's Test P-Value : 0.37109        
##                                          
##             Sensitivity : 0.5556         
##             Specificity : 0.9000         
##          Pos Pred Value : 0.8333         
##          Neg Pred Value : 0.6923         
##              Prevalence : 0.4737         
##          Detection Rate : 0.2632         
##    Detection Prevalence : 0.3158         
##       Balanced Accuracy : 0.7278         
##                                          
##        'Positive' Class : 1              
## 

# 모델 개선형
df_classifier_modify <- naiveBayes(df_train, df_train_labels, laplace = 1)
df_test_pred_modify = predict(df_classifier_modify, df_test)

CrossTable(df_test_pred_modify, df_test_labels,
           prop.chisq = F, prop.c = F, prop.r = FALSE,
           dnn = c('predicted', 'actual'))

##    Cell Contents 
## |-------------------------|
## |                       N | 
## |         N / Table Total | 
## |-------------------------|
## 
## ==================================
##              actual
## predicted        0       1   Total
## ----------------------------------
## 0                9       4      13
##              0.474   0.211        
## ----------------------------------
## 1                1       5       6
##              0.053   0.263        
## ----------------------------------
## Total           10       9      19
## ==================================

confusionMatrix(as.factor(df_test_pred_modify),
                as.factor(df_test_labels), 
                positive = '1')

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 0 1
##          0 9 4
##          1 1 5
##                                          
##                Accuracy : 0.7368         
##                  95% CI : (0.488, 0.9085)
##     No Information Rate : 0.5263         
##     P-Value [Acc > NIR] : 0.05192        
##                                          
##                   Kappa : 0.4633         
##                                          
##  Mcnemar's Test P-Value : 0.37109        
##                                          
##             Sensitivity : 0.5556         
##             Specificity : 0.9000         
##          Pos Pred Value : 0.8333         
##          Neg Pred Value : 0.6923         
##              Prevalence : 0.4737         
##          Detection Rate : 0.2632         
##    Detection Prevalence : 0.3158         
##       Balanced Accuracy : 0.7278         
##                                          
##        'Positive' Class : 1              
## 

3.3 군집화

위에서 구한 유의미한 변수들만 뽑아 군집화를 진행하였다.

library(NbClust)
# 위계적 군집분석 방법
df.scaled <- df_std[,c(17, 4, 31, 5, 36, 22, 24)] # 앞에서 유의미한 결과로 나온 열 선택
#거리 계산
dfDist <- dist(df.scaled, method = "euclidean")

# 다섯가지 옵션 중 이게 가장 예쁘게 나뉘어졌다.
hc <- hclust(dfDist, method = "ward.D2")
plot(hc, hang = -1, cex=.5, main = "Ward Method")

devAskNewPage(ask = T) #새화면 나오기 전에 사용자 대기
nc <- NbClust(df.scaled, distance = "euclidean", min.nc = 2,
              max.nc = 15, method = "ward.D2")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 4 proposed 2 as the best number of clusters 
## * 5 proposed 3 as the best number of clusters 
## * 1 proposed 4 as the best number of clusters 
## * 1 proposed 5 as the best number of clusters 
## * 1 proposed 6 as the best number of clusters 
## * 3 proposed 7 as the best number of clusters 
## * 1 proposed 10 as the best number of clusters 
## * 2 proposed 11 as the best number of clusters 
## * 5 proposed 14 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************

devAskNewPage(ask = F)
par(mfrow = c(1,1))
table(nc$Best.nc[1,]) #각종 추천 알고리즘 결과 요약

## 
##  0  1  2  3  4  5  6  7 10 11 14 
##  2  1  4  5  1  1  1  3  1  2  5

cl.Num <- 2 # 세 개의 군집으로 분류
hc.result = cutree(hc, k=cl.Num)
table(hc.result) # 군집별 데이터 개수 확인

## hc.result
##  1  2 
## 59 14

#군집 결과 시각화
plot(hc, hang=-1, cex=.5, main="Ward Method")
rect.hclust(hc, k=cl.Num) # plot 함수로 시각화 , 지표의 추천 클러스터 개수를 지정

# 군집 결과 시각화2 (using MDS)
mmds2.out <- cmdscale(dfDist)
plot(mmds2.out, type = "n")
text(mmds2.out, rownames(mmds2.out), col = hc.result)

# 군집의 특징 확인
aggregate(df.scaled, by=list(cluster=hc.result), mean)

##   cluster gdp per capita  old ratio    refugees happy score threatened species
## 1       1     -0.4174168 -0.1696201  0.05349294  -0.3450385        0.001928781
## 2       2      1.7591138  0.7148276 -0.22543455   1.4540907       -0.008128433
##   industry employed  unemployed
## 1        0.06917511  0.09155046
## 2       -0.29152368 -0.38581979

결론

다중 선형 회귀 분석 결과 출산율에 영향을 미치는 변수 4개를 찾아냈다. ‘65세 이상 노인 비율’, ‘유아 사망률’, ‘산업 고용률’, ‘지니 계수(소득 불평등 지수)’가 출산율에 영향을 끼쳤다. 노인 비율과 산업 고용률, 지니 계수가 낮을수록 출산율이 높았고 유아사망률이 높을수록 출산율이 높았다. 분석 전 우리는 성차별 격차 지수나 양육 비용과 같은 변수가 출산율에 영향을 끼칠 것이라 예상했지만 기대와 달리 이들은 무의미한 변수라는 것을 알 수 있었다. 결과 중 지니계수를 보면 빈부격차가 적을수록 출산율이 높아지기 때문에 출산율이 높아지려면 소득불평등을 줄이는 것이 근본적인 해결책이다. 정책적 방면으로는 근시안적인 양육비 지원 정책보다 고소득층~저소득층 모두가 살기 좋은 나라를 만드는 정책을 수립한다면, 자연스럽게 한국의 저출산도 해결될 것으로 예측된다. 또한 트위터 분석 결과 우리나라의 사람들은 영어권의 사람보다 출산과 육아에 대해 더 부정적인 반응을 갖고 있었다. 따라서 국가적 차원에서 출산과 육아에 대한 인식을 긍정적으로 바꿀 필요가 있다.

또한 우리는 선진국들이 가입한 다고 알고 있는 OECD국가에 대해 알아보았다. 어느 나라나 가입할 수 있는 명목상의 가입 조건과는 다르게 암묵적인 조건이 있는 듯하다. 대부분의 지표에서 확연한 차이가 보였다. 또한 유의미한 차이가 있는 조건들을 가지고 군집화를 해보았지만 OECD국가들이 군집화 되지는 않았고, 대부분의 유럽나라들이 분리되었다. 또한 분리된 결과를 보니 유럽 국가가 GDP가 높았고, 노인율이 높았고, 행복지수가 높았고, 실업율이 낮은 등의 선진국이라고 불릴만한 지표들이었다.